{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Bayesian test Mrk335(18)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import math\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy.fftpack import fft,ifft\n",
    "import pandas as pd\n",
    "from scipy.optimize import minimize\n",
    "from scipy.optimize import basinhopping\n",
    "from iminuit import Minuit\n",
    "\n",
    "import emcee\n",
    "from pprint import pprint\n",
    "import time\n",
    "from multiprocessing import Pool\n",
    "\n",
    "import random"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 似然函数 p; D = -2 ln p\n",
    "\n",
    "def twi_minus_loglikelihood(log_A,log_f_b,alpha_H,poisson):\n",
    "    alpha_L = 1.0\n",
    "    \n",
    "    perdata18 = pd.read_csv(\"perlist18.csv\")\n",
    "    f = perdata18['f']\n",
    "    per = perdata18['per']\n",
    "            \n",
    "    model = []\n",
    "    f_length = len(f)\n",
    "    for i in range(f_length):\n",
    "        model.append(((f[i]**(-alpha_L))/(1+(f[i]/(10**log_f_b))**(alpha_H-alpha_L)))*(10**log_A)+poisson)\n",
    "     \n",
    "    \n",
    "    length = len(perdata18)\n",
    "    minus_log_p = 0\n",
    "    for i in range(length):\n",
    "        minus_log_p += (per[i]/model[i]+math.log(model[i]))\n",
    "    \n",
    "    \n",
    "    D = 2*minus_log_p\n",
    "    # print (D)\n",
    "    return D"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<hr>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<table>\n",
       "    <tr>\n",
       "        <td title=\"Minimum value of function\">FCN = 4158.28979154928</td>\n",
       "        <td title=\"Total number of call to FCN so far\">TOTAL NCALL = 191</td>\n",
       "        <td title=\"Number of call in last migrad\">NCALLS = 191</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "        <td title=\"Estimated distance to minimum\">EDM = 0.00021316122282715025</td>\n",
       "        <td title=\"Maximum EDM definition of convergence\">GOAL EDM = 1e-05</td>\n",
       "        <td title=\"Error def. Amount of increase in FCN to be defined as 1 standard deviation\">\n",
       "        UP = 1.0</td>\n",
       "    </tr>\n",
       "</table>\n",
       "<table>\n",
       "    <tr>\n",
       "        <td align=\"center\" title=\"Validity of the migrad call\">Valid</td>\n",
       "        <td align=\"center\" title=\"Validity of parameters\">Valid Param</td>\n",
       "        <td align=\"center\" title=\"Is Covariance matrix accurate?\">Accurate Covar</td>\n",
       "        <td align=\"center\" title=\"Positive definiteness of covariance matrix\">PosDef</td>\n",
       "        <td align=\"center\" title=\"Was covariance matrix made posdef by adding diagonal element\">Made PosDef</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "        <td align=\"center\" style=\"background-color:#92CCA6\">True</td>\n",
       "        <td align=\"center\" style=\"background-color:#92CCA6\">True</td>\n",
       "        <td align=\"center\" style=\"background-color:#92CCA6\">True</td>\n",
       "        <td align=\"center\" style=\"background-color:#92CCA6\">True</td>\n",
       "        <td align=\"center\" style=\"background-color:#92CCA6\">False</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "        <td align=\"center\" title=\"Was last hesse call fail?\">Hesse Fail</td>\n",
       "        <td align=\"center\" title=\"Validity of covariance\">HasCov</td>\n",
       "        <td align=\"center\" title=\"Is EDM above goal EDM?\">Above EDM</td>\n",
       "        <td align=\"center\"></td>\n",
       "        <td align=\"center\" title=\"Did last migrad call reach max call limit?\">Reach calllim</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "        <td align=\"center\" style=\"background-color:#92CCA6\">False</td>\n",
       "        <td align=\"center\" style=\"background-color:#92CCA6\">True</td>\n",
       "        <td align=\"center\" style=\"background-color:#92CCA6\">False</td>\n",
       "        <td align=\"center\"></td>\n",
       "        <td align=\"center\" style=\"background-color:#92CCA6\">False</td>\n",
       "    </tr>\n",
       "</table>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<table>\n",
       "    <tr>\n",
       "        <td><a href=\"#\" onclick=\"$('#yGGtVBTMZf').toggle()\">+</a></td>\n",
       "        <td title=\"Variable name\">Name</td>\n",
       "        <td title=\"Value of parameter\">Value</td>\n",
       "        <td title=\"Hesse error\">Hesse Error</td>\n",
       "        <td title=\"Minos lower error\">Minos Error-</td>\n",
       "        <td title=\"Minos upper error\">Minos Error+</td>\n",
       "        <td title=\"Lower limit of the parameter\">Limit-</td>\n",
       "        <td title=\"Upper limit of the parameter\">Limit+</td>\n",
       "        <td title=\"Is the parameter fixed in the fit\">Fixed?</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "        <td>0</td>\n",
       "        <td>log_A</td>\n",
       "        <td>-4.70169</td>\n",
       "        <td>0.50761</td>\n",
       "        <td></td>\n",
       "        <td></td>\n",
       "        <td>-5</td>\n",
       "        <td>-1</td>\n",
       "        <td>No</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "        <td>1</td>\n",
       "        <td>log_f_b</td>\n",
       "        <td>-1.01351</td>\n",
       "        <td>3.30109</td>\n",
       "        <td></td>\n",
       "        <td></td>\n",
       "        <td>-5</td>\n",
       "        <td>-1</td>\n",
       "        <td>No</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "        <td>2</td>\n",
       "        <td>alpha_H</td>\n",
       "        <td>6.39157</td>\n",
       "        <td>4.58253</td>\n",
       "        <td></td>\n",
       "        <td></td>\n",
       "        <td>2</td>\n",
       "        <td>7</td>\n",
       "        <td>No</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "        <td>3</td>\n",
       "        <td>poisson</td>\n",
       "        <td>3.1733</td>\n",
       "        <td>0.103593</td>\n",
       "        <td></td>\n",
       "        <td></td>\n",
       "        <td>0</td>\n",
       "        <td>5</td>\n",
       "        <td>No</td>\n",
       "    </tr>\n",
       "</table>\n",
       "<pre id=\"yGGtVBTMZf\" style=\"display:none;\">\n",
       "<textarea rows=\"14\" cols=\"50\" onclick=\"this.select()\" readonly>\n",
       "\\begin{tabular}{|c|r|r|r|r|r|r|r|c|}\n",
       "\\hline\n",
       " & Name & Value & Hesse Error & Minos Error- & Minos Error+ & Limit- & Limit+ & Fixed?\\\\\n",
       "\\hline\n",
       "0 & $log_{A}$ & -4.70169 & 0.50761 &  &  & -5.0 & -1 & No\\\\\n",
       "\\hline\n",
       "1 & log $f_{b}$ & -1.01351 & 3.30109 &  &  & -5.0 & -1 & No\\\\\n",
       "\\hline\n",
       "2 & $\\alpha_{H}$ & 6.39157 & 4.58253 &  &  & 2.0 & 7 & No\\\\\n",
       "\\hline\n",
       "3 & poisson & 3.1733 & 0.103593 &  &  & 0.0 & 5 & No\\\\\n",
       "\\hline\n",
       "\\end{tabular}\n",
       "</textarea>\n",
       "</pre>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<hr>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "4158.28979154928\n"
     ]
    }
   ],
   "source": [
    "m=Minuit(twi_minus_loglikelihood,log_A=math.log(0.005,10),log_f_b=math.log(1.7E-4,10),alpha_H=3.8,poisson=0.8,\n",
    "         error_log_A=0.1,error_log_f_b=0.1,error_alpha_H=0.01,error_poisson=0.01,\n",
    "         limit_log_A=(-5,-1), limit_log_f_b=(-5,-1),limit_alpha_H=(2.0,7.0),limit_poisson=(0,5),\n",
    "         errordef=1)\n",
    "\n",
    "m.migrad()\n",
    "\n",
    "pprint(m.fval)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAnIAAAHsCAYAAABWqkA7AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzs3XmcXFWd///3p3pJd9JJZ983SEJYwh42WWQTEIggoCDuw4A4I/P9fnVUnBkcZsbvOAv6m4eKaBxw/bIJKgRQwGEVEBJkJySEQMgCSScknXTWXs7vj9NtL+nqruVWnXtvvZ6PRz2q6tatez/dqXS/+5x7zjHnnAAAAJA8mdAFAAAAoDAEOQAAgIQiyAEAACQUQQ4AACChCHIAAAAJRZADAABIKIIcAABAQhHkAAAAEoogBwAAkFDVoQsol7Fjx7qZM2eGLgMAAGBQzz777Ebn3LjB9quYIDdz5kwtWbIkdBkAAACDMrNVuexH1yoAAEBCEeQAAAASiiAHAACQUAQ5AACAhCLIAQAAJBRBDgAAIKEIcgAAAAlFkAMAAEgoghwAAEBCEeQAAAASiiAHAACQUAQ5AACAhCLIAQAAJBRBDgAAIKEIcgAAAAlFkAMAAEgoghwAAEBCEeQAAEi4TZukRx4JXQVCSH2QM7MFZrawubk5dCkAEqKpSVq3zt927gxdDTC4lSsJcpUq9UHOObfIOXdFY2Nj6FIAJMDOndINN0iLFkm33SY99FDoigAgu9QHOQDIR0eHVF8vfe5z0okn+ucAEFcEOQAAgIQiyAEAACQUQQ4AACChCHIAAAAJRZADAAD9WrzYz1GH+CLIAQAQQ3v2hK5Auvde6fnnQ1eBgRDkAACImeZm6V//1d+H5ly0x2trkx58MPrjViqCHAAAMdPa2vs+TXbskJ54wgc6FI8gBwBASpSrlWvHDmnFivKcCwMjyAEAkBL//d/ZX3MuurWDn39e+sUvojlWKeza1T1Io7lZ+v3vw9ZTSgQ5AABSYu3a7K8tXy79+79Hc5729miOUyq/+Y303e/6x2+/Lf3hD2HrKSWCHAAAFWD37tAVlE8cRvyWC0EOAAAgoQhyAACgbHbtiu5YO3ZI69dHd7yB/OY30jPPlOdc+SDIAQCAsti2Tfr+96M73u9/L91wQ3THG8jzz0svvFCec+WDIAcAQApt3x6/VRmKnTuu7/QqLS3+/qGHpPvuK+7YSUWQAwAghV5/3XcHpsm//7v00kt7b3/mmdy7PR9+OLppWOKAIAcAABJh1y6pqam4Yzz6aPmuqysHglxEli2LXxM2AKTVtm2+BaaSptSIk3feCV3BwMxCV1A+BLmIvPdeuhI+AMTZihX+mqg1a0JXUplefTV0BehCkAMAlFVHhw9iW7aErqQ499wj3Xpr6Coqx44doSuIJ4IcAKCsmpr8Op2PPBK6kuK88or02mvhzv+730lLl4Y7f7m9+GLoCuKJIAcAKKuuKST6TiWRNJnAv0H/+Md4zmtWKh0d2V8r12cpjtfeEeQAACiTDRuk1tbQVcTL009Lzc2D79cV5EJ+/+L4xwdBDgCAMvn+96WXX472mNdfL73xRrTHLKff/ja3wRPt7f6+56TCUS731dwsbd4c3fHKpTp0AQCk1auljRsH36+6Wpo3L57N+wByM1AXYSGamrpDTqktXCiddlr+79u0SRoxIvp63n134Ndz+bnaZeFCvxrGtddm3yeOP3sJckAMPPigNGSI1NAw8H4vvSTNnCkNH16WspACW7b4WezHjpVqakJXg6Rbt66wOeS++13pjDOkAw6Ivqb+dAWu731PuvDC3N6T1C5vghwQEyecIM2YMfA+r79enlqQHjfdJO3ZI51yinTMMaGrqWz5tA6lUSHLYrW2SnfeKV18ceHnjWMrWpS4Rg4AUqyjQ5ozp3xdb8iua9L4qqqwdSTJrl1+ihc+v9klskXOzPaV9PeSGp1zF4WuB0B4L78sLV/e/XzyZOnYY8PVA1SKV16RfvlL6atflerrw9QQx9Gk5RKbFjkzu8nMNpjZy322n2Vmy8xshZldLUnOuZXOucvCVAogjpYvl2prpdmzpTFjmDwU6PLjHxc3unP9eulXv8r++n33+ftyrbzQ3l58N3Xf4Hfttf4ShMHEcSnO2AQ5ST+RdFbPDWZWJel6SR+UdKCkj5nZgeUvDUASTJ8uHXKID3MAvFWrCrs+rcuaNQP/YVTu1rA//ckPYgghjgMiYhPknHOPSXqvz+ajJa3obIHbI+lWSeflekwzu8LMlpjZkqampgirBQCgMnStYBHHEFOoNF1zF5sgl8UUSat7PF8jaYqZjTGzH0g63My+lu3NzrmFzrn5zrn548aNK3WtAJBav/iF9I1vFDb1RKVpa/Nd/YW2VD38sPTkk9HWVIyuaZFy6Xos1Pr10uOPR3e8wUaqDnSu1auTFVrjHuT6+6dwzrlNzrkrnXOznHPfLHtVAFBhtm/3oy2jnEk/rdaulW6+ufBVAh591HcfVpIXXpD+53/Kd74tW7K/duONua00ERdxD3JrJE3r8XyqpHWBagGAipb2+biiVskjKcst6u91z2XA4i7uQW6xpDlmto+Z1Uq6RNLdgWsCAKBiDRSa7rzTd012uf9+acmS0teUi4cf7n7c3JyeCZpjE+TM7BZJT0maa2ZrzOwy51ybpC9Iul/SUkm3O+deCVknAAD9aWnxK2ls2xa6ktIaqAv0pZekN9/sfv7UU9LixaWvKRdPP937+UsvDbzurXPSG2+UtqYoxCbIOec+5pyb5Jyrcc5Ndc7d2Ln9Pufcfp3Xw/3ffI9rZgvMbGFzc3P0RQMA0KmlRXr77dyCXEuLtGFD6WtKihATCT/6qHTLLdlfb2qSfv7z8tVTqNgEuVJxzi1yzl3R2NgYuhQAACRJ994rff/7oavobfly6a23ojnW7t3S0qW57x/qmrRs61cn6XrQRC7RBQBAkhUzQW8hBupC7HLzzdLIkdL//t/FnWvJEt8qmU+XaimnNkm71LfIAQDwyCN+TdCBPPmkX8UgSV57Lbf9BupC7CmK+dO2bi38urhiA11Hh58iJ0ktasUiyAEAgnj9denXv/bXi5XaI49Izzwz8D4PPJCONXq//e2957CL6wjNvi2TfQck5KprJO2TT0rXXVdcTX1961vxno6EIAcACGLHDj8R7Ht9F2cskUwOv/GS0JLz298OHHi2bu0/HN91197brrtOeu656GrL1/XX935e7ITTv/99NKGr5+dg27Z4L+mV+mvkzGyBpAWzWUUbAJACTz8tFbLqZH9hpKXFL48VSjlaY9Mu9S1yjFoFAKRFV5dpXFoOd++O9nhPPBHt8XIR59a2XKQ+yAEA4sO5yltHtK/duwuf5qOpKdJSilZM8MplepIbbijs2P2tpZptRYqQLZJRIMgBAMrGOT/oYN680JWEk/Tg0FMx06g8/vjg++zYUdixq6q6Hw/WepnLv0d/1xfGBUEOAFBWZlJSLltua/MtR0nvfhvImjV+gESa9Gx96/m4v5bQ/ubYM+s9YfOrr0ZWWuQIcgAA9LFliw9va9dKt90m3XFH6c/Z2hr9NWe5WLNGevDB8pzrhRe6H+cySfFguo6Rrdu0r5/8pPhzxk3qR60CANKnlAvTt7ZK//Vf0jnnSOPH+23lGFxwxx3hJiQudtqPXD3ySPfjsWOlTZsKP9auXdmnGukZ7OIyMKRUUt8iZ2YLzGxhc3Nz6FIAAHnYskX64Q/37tZqa5Meeqh05+0KAVG0GOVj40Zp+/bynOeHP+z/te3b47NQ/GBTk6S5uzsfqW+Rc84tkrRo/vz5l4euBUBlaGuTFi70LTs1NdLll/t75GfbNumdd/Yeqemc/35OmxamrmJFsQxWMbZs8d/XbN54o3y1DOTJJ8tznv66ZX/zm/KcOwqpb5EDgHJra/O/LD/1KX8Reehf3IiXwdZ8jYM4THPC/5vcEOQAoAQyGWnUqPRfn4P8VSegL6yQZa6Kud4tav21soUYSFIOCfg4AUB6NDf7a3+qqqQJEwh6SI5s1yV2hb6bby5fLVHasyd0BcUhyAGDeO653kPmczVhgvTBD0ZfD5Lt1lv9L7733pO+8AXfagfEwYgRA88n9+yz/W+fNElatao0NXXZuHHg1/u2wPVsfRtsapIHHiisprggyAGDWL1amjxZ2m+/3N+zdasfZk+QQ18dHdJFF0m3317+UZEo3IYNvkXqlFMqoxW1o0Naty63fQf7fnStD1uMH/5QOvnk3PfvGd6eeab488cZQQ7IwZgx0syZue8fp2tFkA6trX5y2kxGmjrV36N8Nm6UHnvMh4lKCHI/+lH342JXfYii6zLXCX/72z/ts4+lPsiZ2QJJC2YnZT0YAOjH0qXS/ff7btlLL5VmzAhdUfk557vBpkwJXUllKWY9VZRe6v+mc84tcs5d0djYGLoUACiYc3590kmTKrdLtrVVeuqpgafvaGryLZeobPm24CVZ6oMcACCcDRv84uO5Xqd0++3Sr3898D4DdW3efnvvbsGBvPZa+haLh0eQAwAM6qmnpHvv9beBZsqvZFu3+jC3YUNu+7/6qrRsWeHny/UXeGurH0H8/PO9ty9dGs3F+VH6wQ+k3/0udBXJUo6lzuIi9dfIAUiG9eulG2/s/Yu4oUH6X/8rXE2Deeopaf586a23pJUrfbcnihey63jr1uzTbBRjxw7/GSnEu+9WxgCLnm65Jb9BEpXUAtcXQQ5ALOzcKU2cKH3yk/55e7v0n/8ZtqZcHHpoemeMr1SlCAUvvlhcQMw3yO3Y4efATKply6T6+miOtW1bNMeJK4IcgCD27Ok9Gm77dv/LqmtxeabXSK7f/tb/8rzootCVxEe5W4zeeqv0k/SG9tvfdj8e6Pub9sFBBLmINLz8RzW++Jw05ijpkEOk2trQJQGxduedfrLlnutOzpkTrh5E57XX/Nxd7e2hKymNlhbp9ddDV1F5+rZKvvxymDrihiAXFdehEcuXSJ/9vvTGG9LBB0tHHSUdfbS/nzuXJgagh7Y26YIL/JQaQJK88YZ0111+vdyBOOdvlXZ9W6lMmBC6gnhKfZAr14TALQe/T+tmvE9Tz5T/c+255/x4+/vuk/7pn/y04Ece6UNdV8CbNo3/4QAQY0uXFv7eH/5Qmj69+7pPFIe2kP6lPsg55xZJWjR//vzLy3bShgbpxBP9rcvGjdKSJdLixdJPf+pXy3auu8Wu6zZ2bNnKBIBKtHu3HwyQi6Ym3xK0fn3+5+laVq0U2tqKP8bixcUfA+GlPsjFxtix0lln+ZvkQ9yaNb7VbvFiPzzv2Wf9op49u2SPOMIHQwCocEuWDPz6T3/q/2YezF13+fnqclFV5UdTFxLkovTuu/7XxHnn+eetrb1ff++9/I+5bl3xdSE8glwoZr5rddo06cIL/baODmn58u5w98tfSi+9JO27b+9wd/DBDKYAUHGWL5dGjMj++ptvdj9etsxfmtyfKBZxL7eOjsqZ5JbVNvJDkIuTTEbaf39/+9Sn/LY9e3yYW7zYB7zrr/ezSvYdTLHfflxAACDVqqqkGTMGHw1bU+NXksgW5IrRtyUsCi0t0qJFyQyYpfDQQ6ErSBaCXNzV1vpBEkceKV15pd/W0iL96U8+3N1zj/SP/yht2uSnmO85mGLqVAZTAKg4+VxqvHy5tGKF9JnP5Lb/N78pffjDBZXVr+Zm39K2bFn3HIqVru+yaXGzY4c0dGjoKroR5JKooUE66SR/67Jxow92ixdLP/mJ9Fd/5UNc38EUY8YEKxtAZbv3Xt/BkO/kuMuWSaNHl6amd9/1AwdCTRrbs3Vx+nQ/tQnibeVKad680FV0I8ilxdix0gc/6G+S/0m5enV3l+x//IcfTDF2bO9wd8QR0rBhYWsHUBG2bpV27cr/fevW+QmkZ82KvqZiZDI+iG3f3vvHaK4DKfo7Xl+33VbYsVA5CHJpZeb/vJs+vfdgimXLusPdbbf5qbFnzdp7MAVt/ABy1Nbmr++qry/uao5ly6KrqZSeesrfT53qB1j0DXKrV0d3rqam3s9vvjm6Y0fBud5L7UWpVMdNG4JcJclkpAMO8LeegylefNGHu6eflr77Xf+T6YAD/HpJs2b1vk2axKAKAL0sWiS98IL06U9L++xT+HFeesm/v7k5utpK4d13/X3f0PrjH3c/HjVK2rzZP371Ven++8tTW7kNNiVMGsXt0nOCXKWrrfWDJObPlz7/eb+tpcW31L3xhr8K+OGHpf/+b/9861Y/HUrfgDdrljRzJi15QAXqGskZxYjOmTP935bZmBXedVlqu3d3P+759y7rsqKUUh/kyrVEV6o0NEjHHutvfbW0+Cs9V6zwwe6ll6Tf/MY/XrdOmjx574A3e7YPf0xsDKBI++/vgxHdbtEaOzb8pMcoTOqDXJAlutKsoUE65BB/62vPHmnVKh/qum5/+IO/X7nSz+TZN+B1PR47Nn7t1QBip7bWr7TQc/JfFK+uLnQFKFTqgxzKqLbWX1c3Z87er3V0SO+80zvkLVrU/bitzQe6ffaRpkzxLXuTJ3c/njLFB0HCXiK99pr0+OO9t23cyOWWSfTCC9K4cYU3sB9zjO8+7dkNWQpRrEUKJAFBDuWRyfgwNmVK7/nvumze7APdW2/5Ltq1a/2FMF2P163z4/z7hru+gW/yZP60jKH166Xx4/281l3M/NgZJMdxx/nr19auLXzVhH328d2jpZz0taFBevRR6fjji1vNsJD1S4FyI8ghHkaN6h50kc22bT7Q9Qx3b70lPflk97Z33vE/xXuGu0mTfIoYN677Nn68785lzdqyGT7cT9eA5Npvv2QstH7MMdKaNYMv5TWYp5/2XzOiwRJkpUGQQ3IMH+6bAQZqCujo8MuV9Q18K1b4yZ+amrpvGzf6yZ96BryeQa/zceMb41SbGSM1j/Q10B8IlFRLi/Tcc+HHRw0fHvb8afP006ErSCeCHNIlk+kOY4ceOvC+HR3Sli29w13XbdUqP0FSU5MOW7FBdTs3S9s3+0X2RozwLYgjR3bf93w8apRqq0ZqxqujpKd6b9eQIeX5PgAJtnu3dNddviu3EFVV0dYDxBlBDpUrk/ELOI4ePWAr36N3+17aI4+Uv4K6udkHwM2b+79fs0a167fosKWbpef7vF5V1W8IPPbdURq5ZLg0dbhvBmho6Pe+vmW4tLNBaqhj4AfQR22t9Bd/IU2YELoSlEKUK2akCUEOyEd1tTRmjL8NoGWTdNf/k/7mb3ps7FrLpp8AuPW3mzVp2Lbu6wBbWvzjrvvOx5/dsE31/9biA2WWsKeGBmno0O5bfX3v57lsGzKEoIjEMfN/HyXpo7t2LQ31ufrDH0JXEE8EOaBczLqD0pQpvV56tVWadJo0asbAh/j+ddLnPicNr2vNGvbU0uK7gHfu9Pc7dvhBID2f93297/PWVj/6t2fQGzKk962uLudtE98aovbqIdLGQd5bW+tXB6mpkapqVdVaI7XX0FeG1Lrzzr1+HOTtoIOkV16Jpp446+gIXUE8EeSAJKqp8U0Po0aV5vjt7dKuXb2D3u7d3bdduwZ/vm2bH1Cye7dGvLVbtnuXtGKQ9+zZ40Nka6sye/boK7tapW+0dn/N2W61tTp1R41qh9VII2o0tqNGZzfXSL/uf99ez6urfVCsqpKqq1XTXqVjllRJrloT1lWprqVKWtV7n677uc9Vqfo3VRq3slpDhlZJW/z2THuVpr9RJT1erSmrqpRZXCU1+PeNXlulquXVGrmhSpk3q6SWHseuqvJd/j1vZlImI9uVUaY1I2vLSO0ZyWWS1fSErKJY2gyViyAHYG9VVX5E77BhkRxu+aM+G048Nff3dLRL3/xX6Zpr5N/cGfD6ve3Zo2d/16rpk1o1Z2artq5r1QtPtGrqGf3v2+t5e7vvqm5v97fdu1W7u13a2K6apjYN3dYuLe3xeo99Z7/Sptqt7Zq8sV3VapMe9dur9rTphFXt0mvtOuXtNtU80y45/94PNLVrxO3tWrClTcN/3r1d7e2+ySHLbV6H00HtHTLnb39unugb/Prc/npPRpnqjKwqIw0ZeN+BbudvzGjUjzI6fVtGw4ZnpGGdQdJM728yzd9hmvCAScNMF7xtam03VT1lmrTL9PF1phH3mtraTY2NJg037bPVdMm7JmemsX/07zt2rWlInWlsk98+5hHT5E0myT9vGG4ae6fphCbTodv983F3mzTC1/HBF03Vr5tUbTrzBVPN6ybV+NfOXmwa0WjSEtP8ZabZzX67ZUwdzpRpNp3xtMnJNPElU+O7ptohpmENpuoa05gNfn8n+3OAHvm8NH2zaViLpGHSyY92B+v6oaYdO3zuH/uiacQ7fruTqaqqc2qUzuMMGya1bDeNeVWauFG9ArpTj7DeY3vtEPvzpMoTVphGru+9/4w1kr3d/3HcetNxK6XJq01D1vhtUzaaOtZIu3dlP2fU2yeskRrWR3/8vnrt11cB7xs+bqR00ILsxywzghyA+OtqrRpgsufm5dLO2ZIOkfasldZulnR2/qdq3S49/n3p2C9La5b4Xul9s/zMvvfb0mWXSS8843ufjz++8xi7pJv/S7r6aukX/y5ddZV/XZJuu0G64ALp9tulSy8d9HLLP3vpBb/SXXOz9P73+4l15dyA4U8dHbrxOx2aMa1DE8d36KgjB9432+3Rhzu0ckWHzlvQoecXd2junA5Nn9Luz++clj7s9PYqp1NPcRo2xenZO512bHeaebFT83qnZx5ymjPbaecOp1n7OjVMcXpvpdPzf3QyOY09xWnoWKfVTzk1Dnd69WW/fc5spxWv+3OYnMaPdxo2y+md5U4bm5zGjXUaMdtpxDi/z7rNTgcd6VRV7bS+yenAg51U41/b+JaTG+2kOU7btjhtyfjtVRmnjnYnN8mpZZg/z54RTruanVytU229k6qd2mq66+jizMmZbxyVdW+XczLnZE7KOKdMu5Tp6Pz3klQlydq7969udapplap3S7WtvY/Tped55Zxq26XMHr+9dps0dEfv/Wu3ODW09H+cuo3SiGanoUOkxubObe9Kje9Je3ZnP2fU20fucKpqLvA4PR6b67m9r+yvFfq+oc9NkUSQA4DUaWmRnniijF1lZt0hV76x8Vd3+daeM8/0c17vGibtaZRaR0sqcCWNPz0mnfNFadQcaeNWafoRknoM9N6wUVpZJR17vKT9pFWv+MDpzpJ2rJJeXyM1zve97RMOl7S/1Py89FrnL/HjzpY0XVorac84qetyr+HHSS/1yO6zZkkTj5PeeMKvtbrvvtLE46Xxs/zrLzRJp/2lpCF+wPhJV0iq96890+TD7zGflpb9wk8tKfkWs7Y26bQvSk/u8NsOOcSvYDF8uJ9PvL7eL03W1777+nB90F9JDeOlR67tfm3MGD+lZV2ddMABfl48yT/v6Og9Oe748dKGDX51w9dfz+3fZPhw//2U+r9G7rjj/NSZ/Tn5ZOmRR6TDD++u65RTpGeflbZuze38UZg7V1q2rHzni8pHPyqNDV1EDwS5iGQy/j/6ypWhK4nGnDnS6aeHrgJIlvXrpeXL/QIl555bnnN2dEiPPebD4z77+CkaRo3ylyeOjfC3zYQJlXNJ3osvhq4AcbZzZ+gKekt9kDOzBZIWzJ49u6TnOfJIvxB0GqxZUxkjoIBSaGwsfCLbQuza5VsBu+ZOq6mJ7NJG5GjjRt+qljRdrXHIz6JFvdeNDi31Qc45t0jSovnz519eyvNUV6dnEsrt20NXACAfNTW+C7CcWlr8ovKTJ5f3vHG0cqV04IGhq0ClSn2QAwBE78EH/eUkF17YvW3bNj84ZMBryFNmxozK6XJGPBHkAAB565r9pGdoW7SosGM9+KA0b17xNQ3mzTfLezE/UA4EOQBAUKtX+1Gc48aV7hz19f4i9ZaWwfcFkiQTugAAAErdHVtVJY0YUdpzACEQ5AAAidTWJi1eHLqKdGA54+QiyAEAEul975Neey10FenAgI3kIsgBQIVbudJPgttztYEkmDat+Dnz2tqiqSWb3buTM6XTrl2hK0AhGOwAABXuV7/yy3rV1Un77Re6msKsWqU/LyIfJ8751T7ibs+e5AV5eLTIAQA0alToCnLX3r73tldflaZPl2prcz9ONU0Zf9bf9xTJQJADACTK8cdLZ57Ze5uZX12hmGu9qqulr3yluNqSgOvh0oW/RwAAJbFnj9TaGv1x587191Ffe2YW71a6qqpoWs4yNOGkCv+cAIBeNm2SrrtOuumm4o7zn/8p3XZbNDVFYetW6cYb/SL3hTrttO7HK1cWX1NPb7wx8OuzZ0d7PqQDQQ4A0MvOnb71Z9OmsHU89VT0x2xvl7Zs8Y9LPWI1Xx0d0qRJoatA0hDkAAB7CX0dVTlan+rqSn+OvsaP33vbhg3dj4cOLey4dJdWLv7pAQB5Wb9e2ry5tOcYObK0x5ekk0/uflxMd2uujjxSGj26NMeurpbOOKM0x0a8xfiyTgAon9tv9wGloyM+yxXdd5+0dKnU0CAde2zoaro99pg0dao0eXLoSqLz0kuhK4jWQIMiHnqofHUMZLBrApEbghwASFq3TvrQh6TGRmnIkNDVeJs2SaeeKt1zT+hK9nbUUdKyZaGriEa2buRt26Rx48pbS1TefDP7a86Vr46BxO0axaQiyAEJYib98pcDT5Ewdar/5Y/8jRoVv4lxC71mqhxqa6W77y7flB1bt5bnPD2tXCkdemj5z1usOE+jgmjxTw0kyCc+MfDcWVu3+l+szzyTfZ8TTvA3oFhnny3t2OFXVSiHpiZpxIj839fRkf/yXVOmSGPHSi+8sPdrdXWsS4r4SH2QM7MFkhbMZgIepMCECYPvc8AB2btOliwJ06qBdMpkyns94fz5/Y/6HMwzzww+MfEbb/h9urpZ998/+wCIv/1b6RvfyL8OoBRSP2rVObfIOXdFY2Nj6FKAshgyxLcY9HerqQldHVB+bW2Dt0K/+6603365dUn23Wf3bv5AQjipb5EDgGI0NfmWzExGOuQQrj1Kg+pqH+76TnFS6MCGF1/0NyCE1LfIAUChpk/3XXnvvis9+KDKblYJAAAgAElEQVQPdUi+2lrp/PN992ka7N7t17VFZeJvSwCJ4Jz0/PO9L1qfMEHaZ5/SnXP8eOncc/3jNWtKd5581NRIt97qQ+bHPx66GowcKa1aFbaG2bOlFSvC1oBwaJEDkAjt7X5E7ubN/rZqlZ8w95VX/K1r/cy0u/hiP3q51CsrIDdz5uT/nqhbz7qW52JetspEixyA2OrokH7+8+7H1dXSBz/on2/dKj3wgA9xkjR8uDRxYpg6y6m6WqqvL93xd++W3n67dMcP6cMfln7969BVROuoo7qnJGpuDlsLwiDIAYilqirpsst6d6X2DDAjRkgXXVT+utJu2TLpqaekgw5K37JVhx7qg39cusmLVVcnnXOO9MgjoStBSAQ5ALE1dWroCirTPvtIJ55YXJAbOtR3h+c7EW++tm/3S5kBlYpr5AAAkfvc56S//MvB96uv9wMGMgX8Nqqr82vjLlmS/3uBtCDIAQDKrq7O359xhh/AkW3h+oFUVUnHHJPfezo6/LJipRaXhelLYfjw0BWgJ7pWAQBlt//+0jXX+DD2zjvlO+977/n7SZPKc5406lophnkV44EWOQBAEOVcp7VLV0tZLusWA0lAkAOACrB1q/TjH0s/+9ngC8gj/rZtG/j1nTvLUwfCI8gBQAXYutVfG7ZuXfqWc9q5U3r22dBVxEulTJANghwAVIwhQ6LvzqypifZ4hWhu9q2M8+aFrgQoP4IcAKBg55wjfeUrhS1V1VNTk/Tuu4W/f+JEqaGhuBrirpApWpB+jFoFgDJYtsyvj9q1nFJaVFX5yX+L1dbmjzNzZvHHSrKOjv63n3MOEx+jfwQ5oATM/MXIt9yS2/5NTYXNo4XkePBBacoU6ZBD/AS42FvoVrU4zP325pv9bx83bu8gt2tXdOeNw9eOwhDkgBIYNcpPctrentv+Rx7pf8kj3U48URo7NnQVyOaAA/rfPnOm9MILZS2l7J54InQFKBRBDigBM2n27NBVoJyWLZNWrw5dReUx861YUaw2MHFi/9uHDSv+2ECpEOQAoEiHHtodJLKFAZTG4YdL998/+LxqQFoR5ACgSOPH+xvKr2vN1lLINvAAiBMGMwMAKtr48b4lte8ce3PnRn+u5ubojxnCSSeFrgBdaJEDIMlPjfGnP/X/2pAh0vHHM7K2FHbt8ist1NaWtnUJ2e23n3Taaf5xz9GbpViPdePG6I8ZwsEHS3feGboKSAQ5AJ3efFNasUI68MC9X3voIemYY+Ixi3/afO97PjyYSX/7t6Gr2duePXQxAnFGkAPwZxMn+iky+nr00fLXUina2qTPfU764Q9DV9K/b3/bB/ja2tCVICpjx6anZRAEOQDAAFpbpS9/Ofo1WhHOkCGhK0CUGOwAAACQUAQ5AAAqSK4rziAZCHIAAFSQSlkmrlK6kAlyAADExBFHlP4cDQ2571tsGAo5pU7PqWTSLJFBzsyGmdlPzexHZvbx0PUAAOLJOWnp0uT8Uj/llNAVRGvMmNAVpF9sgpyZ3WRmG8zs5T7bzzKzZWa2wsyu7tx8gaQ7nHOXS/pQ2YsFACRCW5u0dm3/8yMCaRCbICfpJ5LO6rnBzKokXS/pg5IOlPQxMztQ0lRJqzt347JNAEBW1dXSlCmhq+jta18LXUH6JaUVtlixCXLOucckvddn89GSVjjnVjrn9ki6VdJ5ktbIhzlpgK/BzK4wsyVmtqSpqakUZQNAYrW3S3/4g1+9AeUV+kL8888v/hjHHlv8MVC82AS5LKaou+VN8gFuiqRfSbrQzG6QtCjbm51zC51z851z88eNG1faSgEgYbZtkx5/3C+APnRo6GrKo1JaaQYzfXrxx5g8ufhjoHhxX9mhvyW6nXNuu6TPlrsYAEib+nrphBOkt98OXUl5nHuutHNn6CqQi+pqf40jBhb3ILdG0rQez6dKWheoFgBAws2eHboClEultL7GvWt1saQ5ZraPmdVKukTS3YFrAgAAiIXYtMiZ2S2STpY01szWSPpH59yNZvYFSfdLqpJ0k3PulYBlAkAiNDX569927JAymXQvy3TYYaErAMKJTZBzzn0sy/b7JN1X6HHNbIGkBbNpTwcQgTVrpFtu8Y937vTX8cRRc7NkJl15pfTcc/55WkUxAjOUSZOkd94JXUU6VUrXakx/BEXHObdI0qL58+dfHroWAMm3fbs0caJ0wQVSVVXYJYgGU1cnjR8fugqENHWq/+MD6RX3a+QAIHaqq6Vhw+Id4lBep50mTZjgH48a5e+tv3kXymzkyNAVoNQIcgCQMGa+q/SnP5VWrx58/0q0c6f0xz9GN9nxYCsx1NdLtbW9t114YTTnrkSXXFL8MSqla5UgBwAJM3q09IlP+JbBDRtCVxNfjz4a3fVnoVdiqDSNjaErSA6CHAAkjJk0Y4Y0YkToSlApzj5bev/7Q1eB/qQ+yJnZAjNb2JzmIVsAgFSIa8vfnDnSKaeU95yV0jVarNQHOefcIufcFY200wIA+lFVFboC7+ijpZqa0FXERxwGiyRB6oMcAADZfOlLfq3ZnmbM8FPMIKxCR9zOnevvK6VFL/XzyAEAkM3w4b2fm0mf/ax/3NRU/npQvL6jh9OOFjkAANCvYuahO+ccPy1LKLTIAcAA2tulm2/uf56uU0+V9tmn/DXlq71damvzjyvlhz6SxSzsZ/MDH5B+8Yvc9j3pJOmxx7qfH3WUtHx5aepCN4IcgIK0tUlvvy196lO9tz/5pO+SSkKQu/lm/zVkMv6C97iOGERlmjRJuvRS6VvfClcDU9zEX+qDnJktkLRg9uzZoUsBUieTkaZN672toSFMLYXYs0f65Cel6dPDnH/tWqmjI8y5UT4TJkjr1+f/vtB/XEyeTJBLgtRfI8f0IwDiaNYs6dlnfYistIuzK82cOaEryE3XWrFIltS3yAFAHJ19dvfjXbvC1QEg2VLfIgcAAOLhuONCV5A+BDkAQOTislpCJSl2EuNCruOLs0pZJSPnIGdmdWb2IzM7tpQFAQBK6623pIcf9vdRu+QS6YtfTNagl7Q49NCBX7/44vLUMZByTqWy777lO1dIOQc559wuSZdIqitdOQCAUnvuOemFF/x91OrqohvpeMIJYSeURfQ2bcptvyjWWa2UtVrz7Vp9SNIppSikVMxsgZktbG5uDl0KAMTGqFGhK+jW1QXWd/Tu6af7SWXTbNgwacyY0FWUTybH1FHMihL5Snrrcb5B7npJnzWz68zsVDM7yMwO7HkrRZHFYPoRAMjNihXSH/9Y/vOOHSv93d/5FUEGksbVN4YPl666KnQV6bJ9u7/PtUXufe8rXS3lkO/0I7/rvP9i563nfyvrfM4lrgCQQKtWSePHS+9/f/nPPdhcelOmSPvvX55akGxdgb9SulbzDXKJ6lYFgCR6+23pmWfCjCIcP9539XW1avQnxGoU557rl6zK1Sc/Kd1xR+nqQfmMHi29917oKuIrryDnnHu0VIUAALzVq/3yYWedFW75sGwOOcQPZpg8OXQlA5s1y4fSSlJT41stX3stdCXRmjWLIDeQguaRM7MPmtk1ZrbQzKZ3bjvJzGL+XxsAkmHcOGnevPjNhXXQQdI558RrsAQ8M+mww0JXUXpf/GJu++XTgptkeQU5M5tgZk9LWiTp05IukzS28+XPSrom2vIAAAC6DTYlTde1cSecUPpa4iDfFrnvSmqQtH/nreelhL+XdFpEdQEAgBxUl3HV9JEji19BotQY7DCwsyR92jm3wsz6jk5dI2lKNGUBAIBcHHVU9lUMog4zV16Zffm1PXuiPVcpDR8ubdu29/YDDpCWLi1/PcUo5Bq59izbx0raWUQtAAAgT0OGZB98cvDB0sknR3euurrs120mad79NE1lk2+Qe1zSVX1a47rmkvsL+ZUfYoWVHQAApZatlSq0hgY/aCaUpE3inLR6pfyD3FclHSXpZUn/Ih/iLjezxyQdJ+kfoi2veKzsAAAoFTPpS18afMH6StXWVvh7p0zxo7cxsLyCnHPuZUnzJS2R9Bn5btYLJK2WdIxzbnnUBQIAEGfDh+e+hmi5DbZiRqkVMwXIggXSX/5l/u/L5brAni1vSWyF6ynvsS7OuRWSPlmCWgAAJbJuXe/FwWtq/Nqqo0eHqwmlddllfh3bjRtDV1KY6ur+R+RWVSX764pavvPIXWZmc0pVDAAgevPn+/VTzzije9uBB/qJVa+4IlxdKK1yT7+xeXN5zpPJSJ/6VPbXk97Clq98W+SukzTCzJok/UF+8MPjkp53zgVYfQ8AMJgxY/ytr54tdIifuM/X1ldra+gKKlO+vfqj5a+R+1f5gQ5Xy18vt9nMfmdmfx9xfQAAVJyPftQvhdZ17V3oa92SJJcWuWz7JLE1L9/BDs4595xz7jvOuY845yZJOlPSc5LOkPTPpSgSAJCfDvpI1Noqvfhi6Cpy17MrNJPxz6urpa9/Pfs8ceWsKSneeit0BeWV92AHMztA0ok9blMkvSLpevluVgBAQA0N0gsvhJ0/LA4OO0x69NHQVeTm/POlqVP7fy3kiNghQ/z1aCNHhqsBA8sryJnZBkkjJD0r6TFJfy3pD845ZtsFgJhYsMDfKt2MGT6I7N4dupLBHXaYv4+ya2/IkMH3mThRevfdgffJtvwX4iHfnN8mqUpSbeetpvM5AACIiWuvzW0y3VNOKW6uN4SXV4ucc26ymc2W71I9SX4U60wze02+he5R59xt0ZdZODNbIGnB7NmzQ5cCBGcmvf66dFs//0u3bOEHOhCamXTppb41sVzOP1+64Ybyna+nyZP9+q0rV5b3vEkc1JBNoRMCr5D0Y0kysw9I+ntJV0r6nKRYBTnn3CJJi+bPn3956FqA0ObN87PQZ5O06Q6ANNpvv/Keb8KE8p4P0cr3GrkqSUeoe6DDCfJTkjRLulcMdgBibehQ6YADSnuO+nrpvvv8raczzyzteZPOuXS1EqDyjB4tbdoUuoriJPH/YL4tcs2S6iW9Kz8h8LXy4e0l55L45QOI2imn+BtyV9V5pfE//ZM0d640fXrYeoBCnH++dN115Q9DV10l3XGH9M475T1vXOQb5K6S9Jhz7o1SFAMg+ZI471RoNTXS174mrVolPfRQ6GqQi6OOkhYvDl1FvAwbJo0aJb33XnnPO2ZM9x9DlSjfCYF/3DPEmVlN9CUBABBvJ58cuoLK9NGP9r8931bANPUhFjIh8PskXSN/fdxQM9sh3736L865pyKuDxVk1y5p6dLojzt8uMSgZQBIvgMPjP6YSQ91+Q52+ID8oIZlkv5T0npJEyRdJOkRMzvHOff7yKtERVi5UnrkkWgnn2xt9cf9yleiOyYAoLLMny8tWRK6iv7l2yL3fyXdLekjfQY3/LOZ3SnpXyUR5FCwyZOl886L7njbt0vXXx/d8QAAlWfGjPgGuXxXdjhY0o+yjFBd2Pk6AAC9vP565S1mjjCOO27wfbJ1pyaxmzXfFrktkmZleW125+sAECuvvCJt2+YfZzLS4Yf7kaIoj6lT/bWqpbgGFpB6B7CoR84PGRLvgJdvi9wvJX3TzD5hZnWSZGZ1ZvYJ+W7X26MuEACKdc89UlOTX4bsscek9etDV1RZ6ur8tBTY28iRoStIn1JMgRTnIJdvi9xXJY2R9FNJPzWzFkkNna/d0vk6AMTOaaf5lS3WrAldCeBdc41vIUa0Km0uy7yCnHNup6SPm9m/SDpK0iRJ70ha7Jx7rQT1AQCQSmmcxHbcuPJPCFwOcW6Ry/tvATOrlXSSpPf3uJ3YuR0AACTQ0KHFH+P884s/RrEmTx58nzgHs3zlFeTM7ABJr0u6XtI8Se2d99dLWmFmJZiqrzhmtsDMFjY3N4cuBQCA2DrttOKPUV8vNTQMvl/USjnYwSzewS/fFrmFkpolzXLOHeuc+5Bz7lj5EavNkn4QdYHFcs4tcs5d0djYGLoUAAAiM3Fi6AryN2FC6AoGli2w1ca4zzHfIDdf0tedc2/33Nj5/Ovy180BAFCxynXt26WXluc8UTrllNAVFObAA7Ov8xpavqNW35JUl+W1OklvZ3kNAMomk5E2bpRuvNE/37279+jAJ56Qhg2TNm8OUx/S7UMf8mtH79kTupL4SeqIUjNp9OjQVfQv3yB3taRvmdmbzrmnuzaa2bGS/lnSl6MsDgAKMXGi9KlPSR0d/nlNjZ/LTJJOP93PKSdJkyb5GxCl6dP9/csvh60D2cX5mrd85Rvk/kHSCElPmtkGSRskje+8bZL0d2b2d107O+eOjqpQAMiVmV9NoD8zZ/obAKRBvkHu5c4bACAHdXXSddf5xywLBpRGMS1sPS+7GOg4cW3Fy3dC4M+WqhAASKOrrvK/AMwqexb/rgEAaZwEF8l29NHSAw+ErqJw+bbIAUDFeughqbXVL6Kdq0oObz2deKJ0yCHSmDGhKym/6s7ftJs2FXec+vria4HXs3WtepAkFPcBGgQ5AMjBOed0Lz0U97mw4qi21i/fVIm6Jsjdvr2447zvfdKjj3YP4kFpDRtW/L9ZOfC3IgDkYNo06dBD/S2JE7HmoqUldAUYSE2NNHJkac/BHyndvvjF0BXkhiAHANBRR0lHHinNnRu6EoT0mc+ErqAwUQ1EiGK92XIjyAFAjJn5rrQnnijttTpz5kinniqNH1+6c+Sivl6aMqW71XPGDOmAAyRWWSyPNFyHV0yoO/bY6OooF66RA4AYq6qSLrvMr04xdmzoakpv6FDp8su7n8+a5W9AlLKFvSQOTiLIAUDMZZvcGMjXwQdLs2eHrgJRIsgBAFAhLrzQ33eNwE66ww4r/TniPv1IAhsRAQAAuldLmTOn9OeK66UNBDkAABBbuUwiffrppa+jpsZ3TccNQQ4AAOSttrY85znmGOnss6M95tFH97+9XF9TlAhyAAAgb/vtJ+27b+nPY9bdhRqV6dP73z5mTPd1hElBkAMApFIm4y/qb2sLXUk6ZTL5ze83YoS/P++80tQTlaRNCpz6IGdmC8xsYXNzc+hSAABlNHYsU7fESddkw4cfHr+wtN9+oSsoXOqDnHNukXPuikamBQeAimImjR4duorKc9RR3degDR8etpZc7bNP9tdOOql8dRSCeeQAAEBk3v/+7sfHHis98EC4WorRNX/cgQeGrWMwBDkAObnnnt7L17S3h6sFQDIkccmrpCHIARjURz4ibd++9/Z588pfC+Kpvt63YAwZEroSoLIQ5AAMau7c0BUg7j78Yen882mBQXjOZX/txBOlxx8vXy3lwH85AEiBTEZ67rlwQcqMEJeLuE04W1UVuoLcjB5dvtbegdZWjeMyXbTIAUAKnHqqdOSR3VM8IJ4uvFDasyd0Fd0mTQpdQW5mzJC+9rXQVfgRrMcfH7qK3ghyAJACQ4ZI48eHrqK34cP9dZQzZoSuJD4mTw5dQWXr29q2YEH+76+OWXKKWTkAgLSoq5Muuih0FYWrr/ddj0mZC63SXXCB9NRT+b0n6qW/QiDIIdUyGWnXLumGGwo/xtat2dflA5BeDQ3SNdf0/1p9vXT66dK0aeWtCdkdcoi/vfJK6ErKiyCHVKuvl/76r6XW1uKOM25cNPUASIdMRjrhhNBVpMs110jf+560eXPoSpKFIIfUGzMmdAUAgMEkZQRt3DBYHAAAVIQRI0JXED1a5AAAQEU44gjpgANCVxEtWuQAAMBeBlohIakyGWnYsNBVRIsgBwAA9tLeHrqC8hk6NHQFhaNrFQAA7KWSplaZNy+5Xy8tcgAAoGQuuaS85yukS7iqyq/nmkS0yGEvmYy0dq30ox+V97w7dyZn3T8gyZqaQleASrL//qErKE7fZb3ihiCHvUyfLn3601JHR/nPPXJk+c8JVJIFC/xqJbNnh66kcvC9jq80zF1HkMNeMhkWdgbSKumtI0nU0BC6gsEde6w0Y0boKsrr05+Wpk4NXUXxCHIAAOQoU8SV5XGezuOss0JXUH777BO6gmgQ5AAAyMFf/ZXU2FjYe/fbTxo/Ptp6kqamJnQF6USQAwAgB8UEsUsvja6OJLr00mR0MfcU90EOXQhyAABEIA0XzpfKfvuFriC9CHIAAETgzDOlbdtKe44hQ0p7/HI58UTp/vtDV5EOTAgMAEAEJk0qfcvThz8snXtuac9RDscdF7qC9KBFDgCAhKj0ARPYGy1yAAAACUWQA4CEMpOeey45o+uAJEja/ye6VgEgoU46ya/UMHRo6EqQNF1zulWTAhIvkS1yZravmd1oZneErgUAQhkxwq/jyZJ6yNfw4dLXv56eUbA99VxB48tfDldHuZQ9yJnZTWa2wcxe7rP9LDNbZmYrzOzqgY7hnFvpnLustJUCAJBexSw3lhTDhhV/jI6O4o9RSiH+GX8iqdeqbmZWJel6SR+UdKCkj5nZgWZ2sJnd0+fGmB0AQMUaNsyPXq2rC11JunVdK7drV9g6BlP23nHn3GNmNrPP5qMlrXDOrZQkM7tV0nnOuW9KKnjGHDO7QtIVkjR9+vRCDwMAQGwMG+bXfQWk+FwjN0XS6h7P13Ru65eZjTGzH0g63My+lm0/59xC59x859z8cePGRVctAABADMRlvEp/g31dP9v8C85tknRl6coBAACIv7i0yK2RNK3H86mS1gWqBQCARGI6kcoTlyC3WNIcM9vHzGolXSLp7sA1AQCQKLNmSV/6UugqUE4hph+5RdJTkuaa2Rozu8w51ybpC5Lul7RU0u3OuVfKXRsAAElm5ueIQ+UIMWr1Y1m23yfpvqjPZ2YLJC2YPXt21IcGACAxkrb0FHITl67VknHOLXLOXdHY2Bi6FAAAgpk2bfB9kqKr1XHChLB1xAGXRQIAgEQx80uMVcLqFIPhWwAAABKnXCEu7l3SBDkAAICEIsgBAIDUcFmXE0in1Ac5M1tgZgubm5tDlwIAABCp1Ac5Rq0CAIC0Sn2QAwAASCuCHAAAKMisWVJ9fegqKhtBDgAAFGTePOnjHw9dRWUjyAEAACRU6oMco1YBAEBapT7IMWoVAACkVeqDHAAAQK7iviRXX9WhCwAAoFROOUU64ojQVQClQ5ADAKTWqFH+BqQVXasAAAAJRZADAACp4VzoCsqLIAcAAJBQqQ9yzCMHAADSKvVBjnnkAADAsGGhKyiN1Ac5AACQDKWcw+3zny/dsUMiyAEAkGJDhvj7mpqwdeTi1FNLd+yGhtIdOyTmkQMAIMWGD5euvTZ0FbmZNSt0BclDixwAAEBCEeQAAAASiiAHAACQUAQ5AACALEo5kjYKBDkAAICESn2QY2UHAACQVqkPcqzsAAAA0ir1QQ4AACRDXZ107rndkxgXwrno6kkCJgQGAACxYCbNnx+6imShRQ4AACChCHIAAAAJRZADAABIKIIcAABAQhHkAAAAOg0fLmUSlI4SVCoAAEBpjRkjff3roavIHUEOAAAgoVIf5FiiCwAApFXqgxxLdAEAgLRKfZADAABIK4IcAABAQhHkAABAaji397YkTSeSrxR/aQAAAOlGkAMAAH82e3boCpAPghwAAPizSy8NXQHyQZADAPRrzx6poyN0FSi3NF9PlkbVoQsAAMRPfb106qlSVRW/2IE4I8gBAPaSyUjHHRe6imSrqvL3ZmHrQLrxdxYAACUwd670f/5Pd6ADSoEgBwBACWQyEqtDotQIcgAAINWK6d6Oe9c4QQ4AAJTd4YdLhx4auorkS32QM7MFZrawubk5dCkAAKDTeedJxxwTuorkS32Qc84tcs5d0ciFCgAAIGVSH+QAAADSiiAHAABSw7nQFZQXQQ4AACChCHIAAAAJRZADAAAFGzZMqqsLXUXlIsgBAICCjRolXX116CoqF0EOAAAgoQhyAAAAWWRinpRiXh4AAEAYV14pjR4duoqBEeQAAEAwNTWhK8hu4sTQFQyuOnQBAACgcn3gA9L8+aGrSC6CHAAACGbsWH9DYehaBQAASCiCHAAA6GXaNGnGjNBVIBd0rQIAgF4uuyx0BYVzLnQF5UWLHAAAQEIR5AAAABKKIAcAAJBQBDkAAICESn2QM7MFZrawubk5dCkAAACRSn2Qc84tcs5d0djYGLoUAACASKU+yAEAAKQVQQ4AAKSaWegKSocgBwAxsmOHtGFD6CoAJAVBDgBiYtQov3j4rl3S9OmhqwGQBCzRBQAxMWKEdPHFoasAkCS0yAEAACQUQQ4AACChCHIAACA1xo0LXUF5EeQAAEBqTJ8uXXtt6CrKhyAHAACQUAQ5AACAhCLIAQAAJBRBDgAAIKEIcgAAAAlFkAMAAEgoghwAAEBCEeQAAAASiiAHAACQUAQ5AACAhCLIAQAAJBRBDgAAIKEIcgAAAAlFkAMAAEgoghwAAEBCEeQAAAASiiAHAACQUAQ5AACAhCLIAQCAVDMLXUHpEOQAAAASiiAHAACQUIkMcmZ2vpn9yMzuMrMzQtcDAAAQQtmDnJndZGYbzOzlPtvPMrNlZrbCzK4e6BjOud845y6X9BlJF5ewXAAAgNiqDnDOn0j6nqSfdW0wsypJ10v6gKQ1khab2d2SqiR9s8/7/8I5t6Hz8T90vg8AAKDilD3IOeceM7OZfTYfLWmFc26lJJnZrZLOc859U9K5fY9hZibp3yT91jn3p2znMrMrJF0hSdOnT4+kfgAAgLiIyzVyUySt7vF8Tee2bK6SdLqki8zsymw7OecWOufmO+fmjxs3LppKAQAAYiJE12p/+pvhxWXb2Tn3HUnfKV05AAAA8ReXFrk1kqb1eD5V0rpAtQAAACRCXILcYklzzGwfM6uVdImkuwPXBAAAEGshph+5RdJTkuaa2Rozu8w51ybpC5Lul7RU0u3OuVciOt8CM1vY3NwcxeEAAABiw5zLeilaqsyfP98tWbIkdBkAAKCMFi+Wqqulww8PXUl+zOxZ59z8wfaLy2AHAACAyB11VOgKSisu18gBAAAgTwQ5AGexD3gAAAjLSURBVACAhCLIAQAAJFTqgxyjVgEAQFqlPsg55xY5565obGwMXQoAAECkUh/kAAAA0oogBwAAkFAEOQAAgIQiyAEAACQUQQ4AACChUh/kmH4EAACkVeqDHNOPAACAtEp9kAMAAEgrghwAAEBCEeQAAAASiiAHAACQUAQ5AACAhCLIAQAAJFR16AJKzcwWSFogaauZvd7PLo2SBppkbrDXx0raWHiFsTLY15q08xZ73ELfn8/7ct03l/0G2ofPaXzPy+e0G5/T+J6Xz2m3cn1OZ+S0l3Ouom+SFhb5+pLQX0O5vhdJO2+xxy30/fm8L9d9c9lvoH34nMb3vHxOe73G5zSm5+Vz2uu1WH1O6VqVFhX5epqE+lpLdd5ij1vo+/N5X6775rJfpXxW+ZxG834+p6XF5zSa9/M5HYR1pksUyMyWOOfmh64DGAifUyQBn1MkQdw+p7TIFW9h6AKAHPA5RRLwOUUSxOpzSoscAABAQtEiBwAAkFAEOQAAgIQiyAEAACQUQa6EzOxkM3vczH5gZieHrgfIxsyGmdmzZnZu6FqA/pjZAZ0/S+8ws8+Hrgfoj5mdb2Y/MrO7zOyMcpyTIJeFmd1kZhvM7OU+288ys2VmtsLMrh7kME5Si6Q6SWtKVSsqV0SfU0n6qqTbS1MlKl0Un1Pn3FLn3JWSPiopNlM/ID0i+pz+xjl3uaTPSLq4hOV218eo1f6Z2UnyIexnzrl5nduqJC2X9AH5YLZY0sckVUn6Zp9D/IWkjc65DjObIOnbzrmPl6t+VIaIPqeHyC85Uyf/mb2nPNWjUkTxOXXObTCzD0m6WtL3nHM3l6t+VIaoPqed7/uWpP/nnPtTqetO/VqrhXLOPWZmM/tsPlrSCufcSkkys1slneec+6akgbqkNksaUoo6Udmi+Jya2SmShkk6UNJOM7vPOddR0sJRUaL6eeqcu1vS3WZ2rySCHCIV0c9Tk/Rvkn5bjhAnEeTyNUXS6h7P10g6JtvOZnaBpDMljZT0vdKWBvxZXp9T59zfS5KZfUadrcglrQ7w8v15erKkC+T/KL6vpJUB3fL6nEq6StLpkhrNbLZz7gelLE4iyOXL+tmWtW/aOfcrSb8qXTlAv/L6nP55B+d+En0pQFb5/jx9RNIjpSoGyCLfz+l3JH2ndOXsjcEO+VkjaVqP51MlrQtUC5ANn1MkAZ9TJEHsP6cEufwsljTHzPYxs1pJl0i6O3BNQF98TpEEfE6RBLH/nBLksjCzWyQ9JWmuma0xs8ucc22SviDpfklLJd3unHslZJ2obHxOkQR8TpEESf2cMv0IAABAQtEiBwAAkFAEOQAAgIQiyAEAACQUQQ4AACChCHIAAAAJRZADAABIKIIcgEQxs6+b2Voz6zCzn4SuBwBCYh45AIlhZvPlZ1r/O/l1Nzc4594IWhQABFQdugAAyMP+nffXO+e29reDmdU753aWsSYACIauVQCJ0NmN+vPOp81m5szs5M77M83sbjNrkfS9zv0zZna1ma0ws91mttzMPt3nmGZm15rZBjPbZmY/M7NLO485s3OfrnPM6/PeR8zsjj7bTjCzR81sh5ltMrMfmdnwHq9/pvNYB5vZg2a23cxeM7ML+vl6P2xmz5jZzs5j3WdmM8zsoM5jvL/P/g1m1mJmf1Po9xhA8hDkACTFv0j6RufjUyUdJ2lE5/MbJb0g6UOdjyXpu5L+QdJCSedI+rWkm8zs3B7H/BtJX+/c5yJJOyX9RyHFmdnxkv5H0rudx/rfks6W9ON+dr9ZfuHtD0t6XdKtZja1x7E+KelXkt6Q9FFJn5W0XNK4znUe/9i5raePSKrpPDaACkHXKoBEcM69YWZd18Mtds61mNnJnc9/6Zy7pmtfM5st6fOSPuuc+2nn5t+b2SRJ/yjpHjOrkvRVST90zv1D5z73m9mDkqYUUOK/SXrSOXdxjzrWSvofM5vnnHu5x77/n3Pups59npW0XtK5kn5gZpnOY/3aOfexHu+5u8fjGyX9l5l9wTnX0rnts5IWOec2FlA7gISiRQ5AGtzb5/lpkjok/drMqrtu8i1mh3WGuGmSJkm6q897f5Xvyc1sqHwL4e19zvcHSa2Sjuzzlge6HjjnNknaIKmrRW6upMnqvyWvy62d9x/pPP8sSScM8h4AKUSQA5AG6/s8HyupSlKzfJDquv1EvidikqSJnftu6PPevs9zMarzfN/vc77d8t2d0/rsv6XP8z2S6jofj+m8fyfbyTpb4W5Xd/fqZ+S7dH9XQO0AEoyuVQBp0HcepfcktUk6Xr5lrq8N6v75N77Pa32f7+q8r+2zfbSkrm7MLZ01XCvpvn7Ot66/orPY1Hk/aZD9/lvSE2Y2R9KnJP3MOdeex3kApABBDkAaPSTfQtbonHuwvx3MbLV8K9Z56t2S1XcE6ZrO+wMk/anzvdPku0CXS5JzbruZ/VHSXOfcPxdZ+zJJayV9WtKibDs55540s9ck3SRpunxrI4D/v737d+kqCuM4/n521xY3Fwn6HxpcG90anBqidgdDkAQHa2h0Kb5E5lIg2F+QgoPf3R9TZJtuQZs+Ds8lLopfbIpzfb+We+E+HC53+nCec8+5ZwxykgYnM48jYoP6G3QdGFOty0fAbGY+y8yL7tnbiDgHdoF5KrD1x/oVEQfAakT8oZakLFGzfn2L1I8Nl8AX4DcVsJ4ArzLz5I7vfhkRi8BmRGwCW9Rs3xywlZnjXvl74A2wn5lHd/s6kobENXKShuoltWXJAtXuHFGh6nuv5h2wBjwHvgJTVCC77inwE/jU1b+mZs7+ysw94DHwgNrvbqcb65Sba/gmyszPVKh8SIXCj9392bXS7e764V/GlzQcHtElST3dPnM7wExm/vjPrzNRRLyg9r2bvu2kC0nDZmtVkhrTnToxS7V4R4Y46f6ytSpJ7VkBvgGHwPLkUklDZmtVkiSpUc7ISZIkNcogJ0mS1CiDnCRJUqMMcpIkSY0yyEmSJDXKICdJktSoKzqQKId3jzwOAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 720x576 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "perdata18 = pd.read_csv(\"perlist18.csv\")\n",
    "f = perdata18['f']\n",
    "per = perdata18['per']\n",
    "\n",
    "alpha_L = 1.0\n",
    "log_A,log_f_b,alpha_H,poisson = m.values[0],m.values[1],m.values[2],m.values[3]\n",
    "\n",
    "model = []\n",
    "R_obs = []\n",
    "T_SSE_obs = 0\n",
    "f_length = len(f)\n",
    "for i in range(f_length):\n",
    "    model.append(((f[i]**(-alpha_L))/(1+(f[i]/(10**log_f_b))**(alpha_H-alpha_L)))*(10**log_A)+poisson)\n",
    "    R_obs.append(2*per[i]/model[i])\n",
    "    T_SSE_obs += (((per[i]-model[i])/model[i])**2)\n",
    "    \n",
    "plt.figure(figsize=(10,8))\n",
    "plt.loglog()\n",
    "plt.step(f, per, color=\"b\", alpha=0.5, linewidth=1)\n",
    "plt.plot(f, model, color=\"r\", linewidth=1)\n",
    "plt.xlabel(\"frequency\",fontsize=15)\n",
    "plt.ylabel(\"power\",fontsize=15)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
